from scipy import optimize, special
from numpy import *
from pylab import *

x = arange(0,10,0.01)
for k in arange(0.5,5.5):
  y = special.jv(k,x)
  plot(x,y)
  f = lambda x: -special.jv(k,x)
  x_max = optimize.fminbound(f,0,6)
  plot([x_max], [special.jv(k,x_max)],'ro')

title('Different Bessel functions and their local maxima')
show()